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The fractional Laplacian operator, — ( — A)~ , appears in a wide class of physical systems, including 
Levy flights and stochastic interfaces. In this paper, we provide a discretized version of this operator 
which is well suited to deal with boundary conditions on a finite interval. The implementation of 
boundary conditions is justified by appealing to two physical models, namely hopping particles 
and elastic springs. The eigenvalues and eigenfunctions in a bounded domain are then obtained 
numerically for different boundary conditions. Some analytical results concerning the structure of 
the eigenvalues spectrum are also obtained. 



00 

o 

CD 

e 

I 

> 



O 

o 



> 
in 

(N 

o 
o 



X 



I. INTRODUCTION 

Random walks and the associated diffusion equation 
are at the heart of quantitative descriptions of a large 
number of physical systems P, Despite such ubiq- 
uity, random walk dynamics has limitations, and does 
not apply to cases where collective dynamics, extended 
heterogeneities, and other sources of long-range corre- 
lations lead to so-called anomalous dynamics [3, |3j @- 
To describe these situations, various generalizations of 
Brownian motion have been conceived, generally covered 
under the rubric of fractional dynamics [3| . For example, 
a quite useful model of super- diffusive behavior, in which 
the spread of the distribution grows faster than linearly 
in time, is provided by Levy flights: particles are as- 
sumed to perform random jumps with step lengths taken 
from a distribution that decays as a power law. If the 
variance of the jump length is infinite, the Central Limit 
Theorem does not apply [1, 0, H, H, [13] , and the dynam- 
ics is anomalous. Levy flights, which are dominated by 
rare but etremely large jumps, have proven quite suitable 
in modeling many physical systems, ranging from turbu- 
lent fluids to contaminant transport in fractured rocks, 
from chaotic dynamics to disordered quantum ensembles 
% 5, 11, 12, 13, 14, 15, 16]. 

While the concentration C{x,t) of particles perform- 
ing Brownian motion follows the standard diffusion equa- 
tion, dtC{x,t) = dlC{x,t), the concentration of Levy 
flights satisfies a fractional diffusion equation in which 
the Laplacian operator is replaced by a fractional deriva- 
tive as 



§-Cix,t)^^Cix,t). 



(1) 



order a > [1 

involving a singular kernel of power-law form (see Ap- 
pendix |XT]). For diffusing particles, the index a roughly 



is the Riesz-Feller derivative of fractional 
[isl l , which has an integral representation 



characterizes the degree of fractality of the environment, 
and is in this context restricted to a < 2; for a > 2, the 
correlations decay sufficiently fast for the Central Limit 
Theorem to hold, and Eq. ^ is replaced by the regular 
diffusion equation 

Interestingly, the same Riesz-Feller derivative also ap- 
pears in connection with stochastically growing surfaces 
19, 20]. In this case, the evolution of the height h{x,t) 
of the interface is usually written in Langevin form 



d_ 

dt 



h{x, t) 



d\x 



-h{x,t)-\-'q{x,t), 



(2) 



where r](x, t) represents uncorrelated noise of zero mean, 
and with {r]{x,t)r]{x',t')) = 2T6{x - x')S{t - t'). The 
fractional derivative mimics the effects of a generalized 
elastic restoring force. When a — 2, Eq. ^ describes the 
dynamics of a thermally fluctuating elastic string and 
is also known as the Edwards- Wilkinson equation [2l| . 
However, in many physical systems, such as crack prop- 
agations [l^l and contact lines of a liquid meniscus [23| , 
the restoring forces acting on h(x, t) are long-ranged and 
characterized by a = 1. Other physical systems, such 
as slowly growing films in Molecular Beam Epitaxy, are 
better described by a restoring force that depends on 
curvature, with a = 4 (23 |. 

Better understanding of the properties of the frac- 
tional derivative is thus relevant to many physical sys- 
tems. When the domain over which the operator 
acts is unbounded, the fractional derivative has a simple 
definition in terms of its Fourier transform 



d\x\"^ 



(3) 



More precisely. 



is a pseudo-differential operator. 



whose action on a sufficiently well-behaved function is 
defined through its symbol — I?]". Another form of the 
operator, given in Ref. [1^, is 



d" 



-(-A)t, 



(4) 
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where (—A) is the positive definite operator associated to 
the regular Laplacian, with symbol \q\'^. For this reason. 
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— (— A)'^ is also called the fractional Laplacian. (For 
a = 2 we recover the usual Laplacian [itI. Il8|.) 

Thanks to expression Eqs. ([T]) and ^ on an infi- 
nite or periodic support may be easily solved in the trans- 
formed space. However, whenever boundary conditions 
(BC) break translational invariance, Fourier transforma- 
tion is of limited use, and the long-range spatial correla- 
tions (inherent to the non-local nature of the fractional 
Laplacian operator) make the problem non trivial. 

In this paper we investigate the fractional Laplacian on 
a bounded 1-d domain with various BC on the two sides 
of the interval. In particular, we shall study absorbing 
and free BC: the former naturally arise in the context 
of Levy flights in connection to first-passage problems 
[13 . [2a |. while the latter arise in the context of long- 
ranged elastic interfaces with no constraints at the ends 
[13] ■ The remainder of the paper is organized as follows: 
in Sec. mi we recast Eqs. ([T]) and ^ into the eigenvalue 
problem for the fractional Laplacian. We then introduce 
a specific discretization of the fractional Laplacian, and 
present the main advantages of our choice. In Sec. IIIII 
we discuss the implementation of free and absorbing BC 
by appealing to the examples to Levy fiights and fluc- 
tuating interfaces. The numerical results are presented 
in Sec. IIVI with particular emphasis on the behavior of 
eigenfunctions close to the boundaries. As discussed in 
Sec. |Vl some analytical insights into the problem can 
be achieved by examining certain exactly solvable limits, 
and by perturbing around them. We end with a conclud- 
ing Sec. IVI[ and two short appendices. 



II. MATRIX REPRESENTATION OF THE 
FRACTIONAL LAPLACIAN 

Consider Levy flights in a domain ^ TZ: by apply- 
ing the standard method of separation of variables, the 
concentration C{x,t) in Eq. ([1]) may be written as 

C{x,t) = ^^fe(x)e^'=* / My)C{y,0)dy , (5) 
k 

where ipk{x) and Afc satisfy 

- (-A)t^fc(x) = Afe(a)V'fe(a;), (6) 

with the appropriate BC on dVl. Here — Afc also corre- 
sponds to the inverse of the time constant with which 
the associated cigenfunction V'fc {x) decays in time. Anal- 
ogously, in the context of stochastic interfaces, the shape 
h{x,t) may be decomposed into normal modes h{x,t) = 
J2k hk{t)'4'k{x), where "ipkix) satisfy Eq. ^ and hk{t) are 
time-dependent coefficients. Substituting this expression 
for h{x,t) into Eq. the normal modes are decoupled 
from each other, easing the computation of correlation 
functions. 

For the case of an unbounded domain or periodic BC, 
the set of eigenfunctions and the corresponding spectrum 



of eigenvalues of the operator in Eq. ([6]) is known explic- 
itly By contrast, analytical study of Eq. © with 
different BC is awkward and not completely understood: 
for absorbing BC it has been proven that the operator 

— (—A) 2" on a bounded domain admits a discrete spec- 
trum of eigenfunctions and that the corresponding eigen- 
values are all real and negative and can be ordered so that 

— Ai < — A2 < • • • < — Aoo. However, the exact values of 
the eigenvalues and the corresponding eigenfunctions are 
not known and remain an open question (see e.g. Ref. 

and references therein). It is nonetheless both pos- 
sible and interesting to investigate the properties of the 
fractional Laplacian numerically, and at least two major 
approaches exist for this purpose. 

The first approach consists in implementing the contin- 
uum operator in Eq. ([6]) with a finite differences scheme. 
This is the so-called Griinwald-Letnikov scheme, whose 
construction is directly based on the integral represen- 
tation of the fractional Laplacian operator [2^ [s^, [Ml . 
Considerable insight on the behavior of solutions to the 
fractional diffusion equation on unbounded domains is 
obtained by this method, and it has been shown to be 
highly accurate. However, due to some technical difficul- 
ties, it can not be straightforwardly extended to take into 
account BC [3l,[3l,[35. Another finite element approach 
to discretization of this continuum operator is presented 
in Ref. [35i]. 

The second approach is intrinsically probabilistic in 
nature and consists in replacing continuous Levy flights 
representing ^j^^p with a discrete hops on a lattice: a 
transition probability matrix P; ,„ is constructed, whose 
elements represent the probability of performing a jump 
from position I to m. Analogous to Levy flights, the 
jump probability has a power-law tail which after nor- 
malization reads Pi^m = 1/(2C(« + 1)1^ — m\"~^^), where 
^(.) is the Riemann Zeta function. For this reason, this 
process has also been referred to as a Riemann random 
walk [2^ [3^. The matrix Di^m = Pi,m — h,rm is sup- 
posed to converge to the representation of the contin- 
uum operator when its size goes to infinity. BC can be 
taken into account by properly setting the probabilities 
for jumps leading out of the domain. This approach, how- 
ever, has some shortcomings: first, the convergence of the 
discretized matrix to the continuum operator largely de- 
teriorates as a ^2, i.e. when approaching the regular 
Laplacian [2^, [s^, [sS] . Secondly, it is strictly limited to 
the range a S (0,2], due to its probabilistic underpin- 
nings. 

Our approach is the following: we are interested in rep- 
resenting the action of the operator in terms of a matrix 
A such that the eigenvalues and the eigenvectors of A 
converge to the eigenvalues and eigenfunctions of the op- 
erator when the size M of the matrix goes to infinity. We 
start with the Fourier representation of the discretized 
Laplacian, namely —2(1 — cos(g)) (in line with the sign 
convention in Eq. ([4])), and raise it to the appropriate 
power, — (2(1 — cos((7)))~. The elements of the matrix A, 
representing the fractional Laplacian, are then obtained 
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FIG. 1: Implementing BC in a hopping model: for absorbing 
BC the jump from I to site m' outside the domain leads to 
the death of the particle, while for free BC the jump {I, m') is 
rejected. For both cases, the jump {l,m) within the interval 
is accepted. 



by inverting the Fourier transform, as 



Aq 



_^^^,^l-ra) [2(l„cOs(q))]^ 



(7) 



This is the definition of a Toeplitz symmetrical matrix 
Ai,m [(j)] associated to the generator (the so-called symbol) 
(j){q) = [2(1 — cos(g))] ~ . The generic matrix elements 
depend only on n = |Z— m| and ad hoc algorithms exist for 
calculating the properties of this class of matrices, such as 
its smallest eigenvalue and the determinant [sl, [1^ . 
The integral in Eq. ([7]) may be solved explicitly, to give 



A, 



A{n) 



r(-f +n)r(c 
7rr(l + f + 



^ + 1) . .a 
— — sm(-7r). 
n) 2 



(8) 



In the special cases when a/2 is an integer, A{n) — 
(— I)"~"'+^Cq_s.+„, where Ca,^+n are binomial coeffi- 
cients. We remark that A{n) = for n > a/2, as the 
poles of r(— ^ -I- n) are compensated by the zeros of the 
sin(a7r/2) in Ec^. fS]). The off-diagonal elements Ai^„i^i 
are all positive when < a < 2, but come in different 
signs when a > 2. Thus, for a < 2 the matrix A can be 
normalized and interpreted as transition probabilities for 
a Levy flyer with stability index a. 

While superficially similar, our approach has notable 
advantages compared to Riemann walks. The matrix A 
does not suffer from any deterioration in convergence 
close to a — 2, and can in fact be extended beyond 
the range < a < 2. The relatively simple structure 
of the matrix allows to incorporate BC in a straightfor- 
ward manner. It is also suitable for some analytical treat- 
ments, as we will show in detail in the next Sections. 
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FIG. 2: Implementing BC in a model of elastic springs; Mixed 
BC are imposed by removing all springs connected to sites 
with index m" > M/2 {absorbing BC on the right), and by 
pinning to zero all sites with index m' < —M/2 {free BC 



on the left). 



For the case shown here, 



rpcl. 



^Al^rn{hl 



^Ai,m'hf; -B^ = 0. The interface is free to 



^) 5 -^Z,m' 2 

fluctuate at the right boundary and is constrained to zero at 
the left boundary. 



which do not appear in the case of regular random walks 
need to be introduced, such as between "first passage" 
and "first arrival" times, or between free and reflecting 
BC [H, mj. Therefore, a great amount of ingenuity has 
been employed to solve even apparently simple problems 
such as Levy flights constrained to live on the half-axis 

El- 

The matrix A introduced in the previous Section is 
a priori infinite, thus representing the action of the 
fractional Laplacian operator on an unbounded domain. 
Within our approach, BC can be taken into account by 
modifying the matrix elements related to positions out 
of the considered domain in a suitable manner, as will be 
shown in the following. This modification leads in gen- 
eral to a matrix of finite size Af -I- 1. We will study three 
different kinds of BC: absorbing on both sides, free on 
both sides, and mixed (absorbing on the left and free 
on the right), with reference to two physical models. 
The first concerns hopping particles, the second elastic 
springs: both are well defined for a < 2 and absorbing, 
free and mixed BC are easily implemented. In principle, 
the set of rules by which we will take into account BC 
can be extended to an arbitrary a. 



A. Hopping particles 



III. BOUNDARY CONDITIONS FOR THE 
EIGENVALUE PROBLEM 

Due to the non-locality of fractional Laplacian, it is 
not possible to specify the value of the function ^k{x) 
only locally at the boundaries of a finite domain. Doing 
so leads to erroneous analytical results, in contrast e.g. 
with Monte Carlo simulations [ill, IH, Hi, 0] . This also 
implies that standard techniques such as the method of 
images are not applicable [13, H^l- Subtle distinctions 



Let us consider a particle jumping on a 1-dimensional 
discrete lattice, as shown in Fig. [TJ When the lattice is 
infinite, at each time the particle jumps from position I 
to position m ^ I -\- n [n ^ Q) with a probability Tii^m = 
—A{n)/A{Q). For a < 2 the probability is well defined 
if we set 11/.; = 0, as the elements Ai^m all have the 
same sign. This model is naturally connected to Levy 
flights, since as shown before A represents the discrete 
version of the generator of this stochastic process. Let us 
now discuss how to take into account different BC on an 
interval [-Af/2,M/2]. 
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Absorbing BC are imposed by removing the particle 
whenever a jump takes it to a site m outside the interval. 
In the special case of Brownian particles, BC may be 
assigned locally, since their jumps are of the kind Z — > Z±l 
and they must touch the sites ±M/2 in order to leave the 
interval [3,[l3,[32]- Within our approach, absorbing BC 
are implemented by cutting the infinite matrix 11 into a 
matrix of size (M + 1) x (M + 1), thus setting to all 
the other elements. 

Free BC are implemented as in the Metropolis Monte 
Carlo approach: if the sampled m lies outside the allowed 
interval, then the particle is left at its original location 
I. This means that the element 11;^; is the probability 
to stay at I. From normalization, clearly we must have 
11;^; = 1 — J2i^m^i-^"i- Thesc BC differ from standard 
reflecting BC as implemented e.g. in Refs. [H,!!!], where 
particles abandoning the interval are sent to their mirror 
image with respect to the boundary. Free and reflecting 
BC are identical for Brownian particles, thanks to the 
locality of jumps. 

In the case of mixed BC the particle is removed when- 
ever m < —M/2, and remains at I ior m > M/2. The 
diagonal element of the matrix thus becomes 11;^; = 

1/2 - J2m=l+l'^l,m- 



B. Elastic springs 

Now consider a network of springs connecting the sites 
of a 1-dimensional lattice, as shown in Fig. [2l If the 
spring constant between sites I and m is Ai^m, the asso- 
ciated elastic energy is 



(9) 



where hi is the displacement of site I. The elastic force 
acting on the point {l,hi), is 



Fihi 



SE 
Shi 



J2Ai,mihi-h^). (10) 



Such a model also describes the dynamics interfaces with 
long-range elastic interactions. Let us now discuss how 
to take into account different BC on a bounded interval 
[- M/2, M/2]. 

Absorbing BC are implemented in this case by setting 
hjn = outside the interval [—M/2, M/2], thus cutting 
the infinite matrix A into a matrix of size {M -I- 1) x 
(M -1-1). The diagonal elements are now the same as 
those of the infinite matrix. Physically, this corresponds 
to fluctuating interfaces pinned to a flat state outside a 
domain. 

Free BC are implemented by removing all the springs 
connecting sites inside the interval to sites outside. 
The diagonal elements of the matrix are then Ai^i = 
-~J2i^m^i,rn- Thcsc couditious allow to describe fluc- 
tuating interfaces with no constraints at the ends: in the 
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FIG. 3: Absorbing BC: Convergence of the first eigenvalue 
with M for a — 1.8,2,2.2. Dashed lines are least-square fits 
to straight hnes, and the continuum limit Ai(q) is obtained 
for M-^ -> 0. 



past, these BC have been implemented by using reflect- 
ing BC [13, EE El- We think that our procedure better 
represents the physical situation. 

For mixed BC we set hm = for to < —M/2, and 
cut all the springs connecting I with to > Af/2. The 
diagonal elements of the matrix become Ai^i = A{0)/2 — 



IV. NUMERICAL RESULTS 

In this Section we discuss our numerical results, as ob- 
tained by exploiting the above methods. We will mainly 
focus on the behavior of the first (non-trivial) eigenfunc- 
tion of Eq. ^ , which can be regarded as the dominant 
mode, and of its associated eigenvalue, which represents 
the inverse of the slowest time constant. For simplicity, 
in the following we will assume that fl — [—1,1]. Given 
the matrix A, which now is modified as to incorporate 
the appropriate BC, standard numerical algorithms for 
symmetrical matrices are applied in order to extract the 
spectrum of eigenvalues and eigenvectors. Then, to ob- 
tain the continuum limit, the eigenvalues of A are mul- 
tiplied by a scale factor A X{M/L)", where L = 2 is 
the size of the interval. We remark that, since the first 
eigenvalue for free BC is rigorously zero, we focus on the 
first non-trivial eigenvalue. The eigenvectors of A are 
naturally defined only up to a multiplicative factor, and 
the normalization will be specified later. 

Let us first discuss the finite-size effects: numerical ev- 
idence shows that in the case of absorbing BC the eigen- 
values of A converge to the continuum limit Xk{a) as 
M~^. The finite-size exponent appears to be exactly — 1, 
independent of a, while the overall coefficient increases 
with a. These results are depicted in Fig. [3] for the ffist 
eigenvalue: the continuum limit is obtained by extrap- 
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FIG. 4: Eigenvalues with bsorbing (circles), free (diamonds), 
and mixed (triangles) BC as a function of a. Black squares 
mark the exact values at a = 2 and q = 4 (see Section [VA|) . 



dating the least-square fit of the convergence plot with 
M oo. As opposed to Ref. [2^, our method can be 
extended to any value of a and does not suffer from any 
slowing down in convergence as a — > 2. The extrapolated 
value for a = 2 is A = —2.467 • • • , extremely close to the 
expected value of — 7r^/4. 

Finite-size effects are very similar for mixed BC, while 
for free BC the power law convergence for the first non- 
trivial eigenvalue has an exponent of —2 and the slope 
seems to be approximately constant, independently of a. 

To explore the structure of the eigenvalues of A for 
large M, i.e. in the continuum limit, let us define 



Afc(a) - (-Afc(a))< 



(11) 



In Fig. 0] we plot the behavior of A^ (a) as a function 
of a for absorbing, free, and mixed BC. Note that the 
eigenvalues of the absorbing BC problem exhibit quite 
monotonic behavior and actually seem to lie on a straight 
line: we will come back to this point in Section IV Al 
Moreover, the eigenvalues of free BC seem to be tangent 
to those of absorbing BC close to the point a = 2. 

In Fig. [5] we illustrate the shapes of the ground-state 
eigenfunctions of absorbing BC, corresponding to the 
first eigenvalue, for different values of a. The eigenfunc- 
tions have been normalized such that / ■ip1{x)dx = 1. A 
small and a large value of a have been included to empha- 
size the limiting behavior at the two extremes: for a — > 
the eigenfunction seems to converge to the marker func- 
tion, while for a —^ oo to a 6 function. It can be shown 
that the latter limit is approached so that [s^ 

lim M^)= ^/-y + "\ (l-x^)^. (12) 

Typical eigenfunctions for free and mixed BC are de- 
picted in Fig. El In this case the eigenfunctions have 



FIG. 5: Eigenfunctions with the smallest eigenvalue Ai for 
a = 0.1, 1, 2, 3 and 10 for absorbing BC. The horizontal 
dashed line corresponds to the limiting function for a — > 
(marker function). For comparison, we also show for a = 10 
equation Eq. (|12|) as a dotted line. 





FIG. 6: Eigenfunctions associated with the smallest non- 
trivial eigenvalue for a = 1, 2, 3, for free (left) and mixed 
(right) BC. 



been normalized so that their height ranges respectively 
in [-1,1] and [0,1]. 

An important question is how eigenfunctions behave 
close to the boundaries. As a specific case, we focused 
on the case a = 1, and for absorbing BC, our numeri- 
cal results indicate ipi{x) ~ (1 — ja^j)^^^ as x ±1 (see 
Fig. [7|). This result is consistent with the findings of 
Refs. [H, H^, which show that in general for absorb- 
ing BC the eigenfunctions scale as {—\x\ + 1)"/^. The 
limiting behavior for free BC in Fig. [7] is less clear: the 
convergence is rather poor, and we are unable fully char- 
acterize the dependence of the slope on a. Nonetheless, 
we can exclude the simplest ansatz that the eigenfunction 
for a generic a scales linearly close to the boundaries, as 
suggested by the behavior at a = 2 and a = 0, where 
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FIG. 7: Scaling of the first eigenfunction close to the boundary 
for fractional Laplacian of a = 1, with absorbing (top) and 
free (bottom) BC. Symbols correspond to numerical eigenvec- 
tors for M = 256, 512, 1024, while solid lines correspond to 
[x + 1)^/^ and {x + 1)^/^, respectively. 



tpi{x) ~ (1 — l^l)-'^. In fact, the fit in Fig. [7] is for an 
exponent a/2 + 1 = 3/2. 



ANALYTICAL RESULTS FOR ABSORBING 
BC 



For the case of absorbing BC it is possible to derive 
further information on the structure of the eigenvalues of 
Eq. ([6]) by resorting to analytical treatment. 



A. Even a, and general structure of the eigenvalues 

When a is an even integer, the eigenvalue- 
eigenfunction Eq. ([5]) may be cast in a different way. 
In particular, Eq. ([3]) can be extended to complex q by 
omitting the absolute value. Then, since A = — is 
real and negative, we can associate to each A^, a in- 
dependent solutions characterized by qj — AkUjj, for 
j — 0, I,-- - ,a~l, where ujj = cos{2TTj/a) + ism(2TTj/a) 
are the a roots of unity. The general form of an eigen- 
function is 



3=0 



(13) 



where Cj^k are to be determined by imposing the BC 

M±i) = = V-f/'^'^lil) = 0. (14) 



Thus, determining is equivalent to finding the zeros 
of the determinant of the a x a matrix B 



B 



— lAojo 



— iAuJct-1 



, a/2-1 «Aw„_i 
a/2-1 



\ 



J 



(15) 



The structure of the function det(i3) = is rather 
involved. However, for large k it is possible to rewrite 
this equation in the following form 



/a(Afe)cos(2Afc) -|-5„(Afc) = 0, 



when a/ 2 is even and 



/a(Afc) sin(2Afe) + ga{Ak) = 0, 



(16) 



(17) 



when a/2 is odd. Here fai-^k) = cosh(2cot(7r/a)Afe) and 



ga(Afc) 
/a(Afe) 



~2sin(^)Afc 



(18) 



when fc — > oo. 

Two special cases need to be considered separately: for 
a = 2 we have (72 (A^) = and for a = 6 an acciddental 
factorization gives geiAk) — sin(Afe) (cosh(A/3Ai;) -f • • • ) . 
This allows to conclude that for large k the roots of 
det(i?) = converge exponentially fast to those of 
cos(2Afe) = when a/2 is even or sin(2Ai,) = when 
a/2 is odd. These asymptotic roots are exact for a = 2 
for every k and for a = 6 for all odd k, thanks to the 
factorization. 

These considerations, together with the fact that 
Ak{a) < Ak{a + 2), allow to state that the eigenval- 
ues Ak{a) as a function of k will be better and better 
described by a monotonically increasing function whose 
simplest form is the straight line 



AlPP^-{a)^^a + ^{2k~l). 



(19) 



Equation (|19p is consistent with our numerical findings 
and generalizes an observation by Rayleigh, that for a — 
4 the two values Ak{a) and A^^^^'{a) are identical to 
the sixth decimal digit for fc > 4 ;48j|. In particular we 
remark that direct numerical evaluation of det(i?) = 
reveals that Eq. is a very good approximation even 
for A; = 1 if a is not too large, while it has been shown 
that for very larg e a the asymptotic behavior of the first 
eigenvalue is [33] 



Ai{a) = {Aan) 



1 a 



(20) 



Surprisingly, the asymptotic form of Eq. (|T9|) is valid also 
for a generic real a, as shown in Fig. [5] for k — 1 and 
k — 2. Setting aside some special cases of a such as 
2 and 4, to our best knowledge the approximation in 
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FIG. 8: Afc as a function of a for fc = 1 and 2 (dots), compared 
to the approximation in Eq. (I19|) (straight Hnes). 
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FIG. 9: The difference between Afc(a) to A^PP=' (a) for a = 1 
(squares), a = 2.5 (diamonds) and a = 4 (dots), as a function 
of Afe. 



Eq. (|T9|) is a new result. To illustrate the trends, the 
error in the approximation in depicted in Fig. [S) In all 
cases considered, numerical results indicate that the error 
vanishes exponentially for large fc, in agreement with the 
analytical findings for even a. 



B. Perturbation theory 

We next examine the behavior of eigenvalues 
close to a = 2 and a = using standard per- 
turbation theory. Throughout this Section we will 
consider a symmetric domain Q, — [— L/2,L/2]. 



1. Perturbation around a = 2 

The ground state eigenvector for a = 2 on the discrete 
interval [-Af/2,M/2] is 



with a corresponding eigenvalue of 
'M 



Ai = 



L 



(V'll^lV'i), 



(21) 



(22) 



where L is the length of the interval. In order to deal with 
dimensionless quantities, we multiply Ai by L", and set 



Ai = AiL" = M"(^i|A|V^i). 



(23) 



For a = 2, where A{0) = -2, A{1) = 1 and A{n > 1) 
0, we have 



Ai 



2-2cos(-) 



(24) 



Setting a — 2 + e, the operator A{n) becomes, at the first 
order in e: 



A{n) 



-2-e 

l + |e 
1^ 



for n = 
for n ~ 1 
re for n > 1 



(»i+l)Ti(n-l) ' 

The correction to the ground state is given by 

A* = Ai + (5A = Af2+-(V'i|^|V'i), 
which can be rewritten in the following way: 



(25) 



(26) 



A 



M 



M/2~ 



^ = A(0) + 2^^(n) Ml)Ml + n). 

l=-M/2 



By noticing that 

M/2-n 
l=-M/2 



M 



M 



fmT\ 1 fn'K\ 
■ cos — H — sm — , 



we can rewrite the previous expression as 



A* = -A^2^ 



i- 



where Q, in the limit of large Af, is given by 

M 

2 4 A/2 



2 
M2 



"25:^(n)(l--^ 

(ttx) 



n=2 



dx 



(1 - x) cos(^x) + 2Hli££2 - 1 + ^ 



Performing the integration, we find 
QAf2 = _7r2 log(Ar) + TT (Si(7r) + 7rlog(7r) - 7rCi(7r)) 
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FIG. 10: The error in slope of 5Xk, compared to Eq. pop for 
Of = 2 as a function of fcvr (asterisks). The enveloping dashed 
curves are ±4/(fe7r)^. 



where Si and Ci are the Integral Sine and Integral Cosine 
functions, respectively. We can finally come back to A|, 
which, expanding for small e, reads 

A* = -TT^ + e [7r2Ci(7r) - 7rSi(7r) - log(7r)] . (27) 

This approach can be extended to eigenfunctions "iAfcCO 
of every order k. By replacing "01(0 i'^to Eq. (|26p with 
the generic ipk{l) (see Appendix I A 2[) and performing the 
summations as shown above, after some algebra we find 
the first-order correction SXk = A^, — Afe, with 

^Afc = e [fcV^Ci(fc7r) - fc7rSi(fc7r) - fcV^ log(fc7r)] . (28) 

Now, consider the curve A^'^^'^', which after rescaling 
by a factor L" gives 



A appx. 



Ja+|(2fc-l) 



(29) 



By putting a — > 2 + e and expanding for small e, we get 



^^appx. 



2 



fc TT" log(/c7r) 



(30) 



We can thus compare Eq. ((28)l . which derives from the 
perturbative calculations, with Eq. (|5D|) . which stems 
from our generic approximation to the eigenvalues of 



Eq. (O. In Fig. [TO] we plot the error 5Xk 



SXl 



a function of kn. As k increases, the slope of the curve 
along which the actual eigenvalues lie in the proximity of 
a = 2 approaches very rapidly to the slope of the curve 

\ appx. 

We have also applied perturbation theory for a = 2 
to the case of free BC, for which the eigenfunctions are 
known analytically (see I A 2|) . Calculations analogous to 
those leading to Eq. (pS)) allow to derive SXk as 

5Xk = e [4 + k^Tr^Ci{kTT)+ 

-3fc7rSi(fc7r) - Pt:^ log{kn) + 2fc7rSi(2/c7r)] . (31) 



The values of SXk for free BC are close but not 
equal to those of absorbing BC, thus ruling out 
the hypothesis that the curves Afe(a) for free and 
absorbing BC are tangent near the point a = 2. 



2. Perturbation around a = 

When a is 0, becomes the identity operator — / 

and the associated first (and only) eigenvalue is Ai(a) — 
1. In principle, for a = the operator is highly degener- 
ate, but considering the limiting behavior and the scaling 
behavior near the boundaries we are led to conclude that 
the discrete ground-state eigenvector for a = is 



Ml) 



1 



(32) 



where In{l) is the marker function of the domain = 
[-M/2, M/2] (see Fig.E]). Setting a = 0+e, the operator 
A{n) is corrected at the first order as 



A{n) 



-l + o{e^) for n = 



for n > ' 

The correction to the ground state is given by 

T 



XI = 



M - 



-Y,In{l)A{n)In(m), 



which in the limit of large M is 

AI = -Ar [l-elog(Af) + e(l-7)], 



(33) 



(34) 



(35) 



where 7 = 0.57721566 • • • is the Euler-Mascheroni con- 
stant. Expanding for small e, we finally get 



a: 



-l-e(l-7). 



(36) 



This value is to be compared with X^^^- ^ which for a 
+ e reads 



Af 



-l-elog(- 



(37) 



C. First passage time distribution 

Knowledge of the fractional Laplacian operator allows 
us to address the temporal behavior of the Levy flyer 
concentration C(x, t|a;o)i where xq is the starting position 
of walkers at t = 0. For example, let us consider the 
first passage time distribution for the one-dimensional 
bounded domain SI with absorbing BC on both sides, 
which is obtained as [43]: 



d f 

p{t\xo)^-—j AxC{x,t\xo). 



(38) 
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FIG. 11: MFPT as a function of the starting point xo for 
a = 1, 1.5 and 2. Here L = 2 and M = 1024. Solid lines 
are the analytical result {t^){xo) = (1 - Xo)°'^^/r{a + 1), 
while dashed lines are obtained from the numerical solution 
{t^){xo) = -A"^l(2/M)°. In the limit of large M, the two 
results are in complete agreement for all xq and a. 



Analogous calculations for the second moment m = 2 
lead to 

^l^{t^)i:,o) = ~2(t')ixo). (42) 
d\xo\" 

More generally, the moments of the first passage time 
distribution are obtained recursively from 

^^(f")(xo) = -m(f"-i)(xo), (43) 

for m = 1, 2, • • • . 

This above expression can be rewritten as 

Solving numerically this relation, namely {t"^){xo) — 
(-l)"r(m+ 1)(L/M)™"A-'"1, allows us to compute all 
the moments of the first passage times distribution, which 
is akin to knowing the full distribution. 

VI. CONCLUSIONS 



In particular, moments of the distribution p{t\xQ) are 
given by 



{nM= / dtt"^p{t\xo) = 

Jo 



dt 



dtt™- / Cix,t\xo). (39) 



For m = 1, integrating by parts a using the relation 

^^C{x,t\xo) ^ ^^C{x,t\xo), (40) 

we get 



(^^X^o) 



dx C(x,oo|a;o) - / da; C(x, 0|xo) = -1. (41) 
n 



This equation for the mean first passage time (MFPT) 
may be solved analytically in closed form (see Ref. [3^, 
and references therein), to give {t^){xo) = ((L/2)^ — 
Xq)"^"^ /T{a + 1), where L is the length of the bounded 
interval (we have assumed that the interval is symmetric 
around the origin a; = 0). In Fig. [TT] we compare this 
expression with the numerical solution obtained by re- 
placing the fractional Laplacian with the discrete opera- 
tor A, namely (t^)(xo) = —A'~^\{L/M)°'; the two curves 
are in excellent agreement for all a and xq. We remark 
that the required inversion of the discrete operator may 
be efficiently performed thanks to the fact that A is a 
Toeplitz matrix [40|] . 



In this paper, we have studied the eigenvalue- 
eigenfunction problem for the fractional Laplacian of or- 
der a with absorbing and free BC on a bounded do- 
main. This problem has applications to many physical 
systems, including Levy flights and stochastic interfaces. 
We have proposed a discretized version of the operator 
whose properties are better suited to bounded domains. 
It does not suffer from any slowing down in convergence 
and can easily take into account BC. When a < 2, the 
discrete fractional Laplacian may be interpreted in the 
light of two physical models for hopping particles and 
for elastic springs, where the BC emerge naturally and 
are easily implemented. An analytical continuation for 
a > 2 is also discussed. Our approach easily allows to 
obtain the numerical eigenfunctions and eigenvalues for 
the fractional operator: eigenfunctions corresponding to 
absorbing BC show the expected power-law behavior at 
the boundaries. We also gain analytical insights into the 
problem by calculating perturbative corrections for the 
eigenvalues around a = and 2. Further information on 
the eigenvalue structure is obtained by studying the case 
of even a, where a semi- analytical treatment is possible: 
for every a the spectra seem to approach exponentially 
fast a simple functional form. This conjecture has been 
proven for the case of even a and is supported by numer- 
ical investigations for real a. The first passage problem 
and its connection to the fractional Laplacian operator 
were also explored. 
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The associated eigenvalues are = ((/c — l)7r/2)^, where 
fc = 1, 2, • • • . For mixed BC, namely 1) — V'i^Hl) = 
0, we have 



APPENDIX A: ADDITIONAL NOTES 



1. Integral representation of Riesz derivatives 

Riesz fractional derivatives are defined as a linear com- 
bination of left and right Ricmann-Liouville derivatives 
of fractional order, namely 



If /(2fc-l 

71 [-^ 



— 1)ttx 



M-ir^sini^^^^)]. (A6) 



where 



1 



2 cos((m — a)7r/2) 



VI = 



and 



1 



r(a) J a 



(Al) 



dy [x-yr-'^-' f(^\y) (A2) 



and the associated eigenvalues are Afe — {{2k — l)7r/4)^, 
where k = 1, 2, • • • . 

For absorbing BC, we present here also the analyti- 
cal expressions for the eigenfunctions corresponding to 
the first even values of a. For a = 4, the condition 
det(i?) — becomes cos(2Afc) cosh(2Afc) — 1, whose first 
roots are Ai — 2.36502 • • • , A2 = 3.9266 • • • , and so on. 
Correspondingly, the normalized eigenfunctions are 



■D- = TTr-J (2/ - x)™— 1 (A3) 

with a G (to — 1,to), to integer, and x €z — [a,b]. 
This definition does not hold for odd a. The integrals in 
Eq. (jAl|) have a power-law decaying kernel [13, ■ 



2. Eigenfunctions of 



-A) 2 for even a 



When a = 2 the operator in Eq. © is the regular 
Laplacian. For the case of absorbing BC we impose 
V'fe(-l) = -0^(1) = and get 



ipk{x) 



cos(^l^) when k is odd 
sin(^^) when k is even 



(A4) 



The associated eigenvalues are Afe — (A:7r/2)^, where k — 
1, 2, • • • . For the case of free BC we impose = 
V'i^'(l) = and get 



i'k{x) 



sm 



2 

{k—l)'TTX 



when k is even 



(A5) 



ipk{x) 



cos(Afcx) 

y2cos(Afc) 
sin(Afc£)_ 



cosh(Afca:) 



when k is odd 



V2cos(Afc) V2cosh(Afc) 



V2cosh(Afc) 

sinh(Afc£)_ ^^Q^ ^ even 



(A7) 



For the case a — 6, due to a highly symmetric struc- 
ture of the determinant equation, eigenfunctions may be 
expressed in close form. For example, the normalized 
ground state eigenfunction is 



■01 (a;) = tanh ( — ^ ] cos(7rx) 



V3 



cos h(V37r/2) 

1 /TT 

— sm — 

cosh(V37r/2) V2 



cos (^'^xj cosh ( X 



X 1 sinh f -^-^—x 



(A8) 
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